library(logr)
log_open("log_main.log")
log_code()
####################
#REPLICATION FILES: MAIN FILE
#Article: "When Does Online Public Diplomacy Succeed? Evidence from China's ‘Wolf Warrior’ Diplomats"
#Authors: Daniel Mattingly and James Sundquist
#This Version: June 18, 2022


####################
#Description: Running this script will create the figures and tables in the main text of the article.
#The script was written with R version 4.1.2
#It was created and tested on Mac OS X (11.6.5).
####################

##Clean the workspace
rm(list=ls())


#Manually set the working directory here:

setwd("~/Dropbox/Research/Diplomacy/Replication")

#For example, if the replication directory has been unzipped in a Mac Downloads folder, 
#in a folder you have renamed Mattingly_AJPS_Replication
#the above might read something like: setwd('~/Downloads/Mattingly_AJPS_Replication') 


####################

#Check to see if necessary packages are installed.
#If not, the packages will be installed automatically.

packages <- c("tidyverse", "rio", "coefplot")
new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
rm(packages, new.packages)

# Load required packages
library(tidyverse)
library(rio)
library(coefplot)

#Load data
pooled <- import("twitter_diplomacy_data.csv")

# Ensure levels of treatment variable appear in proper order
pooled$Treatment <- factor(pooled$Treatment, levels = c("control", "prochina", "antius"))
pooled$t_ <- pooled$Treatment

#Create data subsets for pre- and post-Galwan Valley Incident
may <- pooled %>% filter(galwan == 0)
june <- pooled %>% filter(galwan == 1)

####
# Main Analysis
####

# Fit models
m.gov.may <- lm(china_gov_pca ~ t_, data = may)
m.gov.june <- lm(china_gov_pca ~ t_, data = june)
m.gov.pooled <- lm(china_gov_pca ~ t_, data = pooled)

m.people.may <- lm(china_people_pca ~ t_, data = may)
m.people.june <- lm(china_people_pca ~ t_, data = june)
m.people.pooled <- lm(china_people_pca ~ t_, data = pooled)

m.policy.may <- lm(china_policy_pca ~ t_, data = may)
m.policy.june <- lm(china_policy_pca ~ t_, data = june)
m.policy.pooled <- lm(china_policy_pca ~ t_, data = pooled)

m.covid.may <- lm(china_covid_pca ~ t_, data = may)
m.covid.june <- lm(china_covid_pca ~ t_, data = june)
m.covid.pooled <- lm(china_covid_pca ~ t_, data = pooled)

us.gov.may <- lm(usa_gov_pca ~ t_, data = may)
us.gov.june <- lm(usa_gov_pca ~ t_, data = june)
us.gov.pooled <- lm(usa_gov_pca ~ t_, data = pooled)

us.people.may <- lm(usa_people_pca ~ t_, data = may)
us.people.june <- lm(usa_people_pca ~ t_, data = june)
us.people.pooled <- lm(usa_people_pca ~ t_, data = pooled)

us.policy.may <- lm(usa_policy_pca ~ t_, data = may)
us.policy.june <- lm(usa_policy_pca ~ t_, data = june)
us.policy.pooled <- lm(usa_policy_pca ~ t_, data = pooled)

us.covid.may <- lm(usa_covid_pca ~ t_, data = may)
us.covid.june <- lm(usa_covid_pca ~ t_, data = june)
us.covid.pooled <- lm(usa_covid_pca ~ t_, data = pooled)

# Coefficient Plots
coefmaker_21 <- function(models, outcomes, term, title){
  outcome <- myoutcomes
  estimate <- unlist(lapply(models, function(x)coef(x)[[term]]))
  se <- unlist(lapply(models, function(x)coef(summary(x))[term, 2]))
  lower <- estimate - 1.96*se
  upper <- estimate + 1.96*se
  df <- cbind(outcome, estimate, se, lower, upper) %>% as.data.frame()
  df$outcome <- factor(df$outcome, levels = rev(outcomes))
  df$estimate <- df$estimate %>% as.character() %>% as.numeric()
  df$lower <- df$lower %>% as.character() %>% as.numeric()
  df$upper <- df$upper %>% as.character() %>% as.numeric()
  
  fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome))
  fig1 <- fig1 + geom_point()
  fig1 <- fig1 + geom_errorbarh(height = 0)
  fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
  fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                            panel.grid.minor = element_blank(),
                                                            panel.background = element_blank(),
                                                            axis.title.y = element_blank(),
                                                            legend.title = element_blank())
  fig1 + ggtitle(as.character(title))
}
# input 1: list of `mymodels`
# input 2: list of `outcomes` 
# input 3: which `term` of the regression you want
# input 4: `title`


# Figure 1 in article
mymodels <- list(m.gov.pooled, m.people.pooled, m.policy.pooled, m.covid.pooled)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
pdf("Figure1.pdf", width=3.5, height = 3.5)
coefmaker_21(models = mymodels, outcomes = myoutcomes, term = 2, title = " ")
dev.off()

# Figure 2 in article
term <- 3
mymodels <- list(m.gov.may, m.gov.june,
                 m.people.may, m.people.june,
                 m.policy.may, m.policy.june,
                 m.covid.may, m.covid.june,
                 us.gov.may, us.gov.june,
                 us.people.may, us.people.june,
                 us.policy.may, us.policy.june,
                 us.covid.may, us.covid.june)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
wave <- rep(c("normal", "crisis"), length(myoutcomes))
estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
lower <- estimate - 1.96*se
upper <- estimate + 1.96*se
df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
df$wave <- factor(df$wave, levels = c("normal", "crisis"))
df$estimate <- df$estimate %>% as.character() %>% as.numeric()
df$lower <- df$lower %>% as.character() %>% as.numeric()
df$upper <- df$upper %>% as.character() %>% as.numeric()
df$country <- c(rep("China", 8), rep("USA", 8))

fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          panel.background = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          legend.title = element_blank())
fig1 + facet_wrap(~ country)
pdf("Figure2.pdf", width = 4.5, height = 3.5)
fig1 + facet_wrap(~ country)
dev.off()

# Figure 3 in article
coefmaker_3 <- function(models, outcomes, term, title){
  outcome <- unlist(lapply(outcomes, function(x)rep(x, 2)))
  wave <- rep(c("normal", "crisis"), length(outcomes))
  estimate <- unlist(lapply(models, function(x)coef(x)[[term]]))
  se <- unlist(lapply(models, function(x)coef(summary(x))[term, 2]))
  lower <- estimate - 1.96*se
  upper <- estimate + 1.96*se
  df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
  df$outcome <- factor(df$outcome, levels = rev(outcomes))
  df$wave <- factor(df$wave, levels = c("normal", "crisis"))
  df$estimate <- df$estimate %>% as.character() %>% as.numeric()
  df$lower <- df$lower %>% as.character() %>% as.numeric()
  df$upper <- df$upper %>% as.character() %>% as.numeric()
  
  fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
  fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
  fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
  fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
  fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
  fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
  fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                            panel.grid.minor = element_blank(),
                                                            panel.background = element_blank(),
                                                            axis.title.y = element_blank(),
                                                            legend.title = element_blank())
  fig1 + ggtitle(as.character(title))
}
mymodels <- list(m.gov.may, m.gov.june,
                 m.people.may, m.people.june,
                 m.policy.may, m.policy.june,
                 m.covid.may, m.covid.june)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
pdf("Figure3.pdf", width = 3.5, height = 3.5)
coefmaker_3(mymodels, myoutcomes, 2, "")
dev.off()



log_close()
